Reaction rate calculation by parallel path swapping 
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The efficiency of path samphng simulations can be improved considerably using the approach 
of path swapping. For this purpose, we have devised a new algorithmic procedure based on the 
transition interface sampling technique. In the same spirit of parallel tempering, paths between 
different ensembles are swapped, but the role of temperature is here played by the interface position. 
We have tested the method on the denaturation transition of DNA using the Peyrard-Bishop-Dauxois 
model. We find that the new algorithm gives a reduction of the computational cost by a factor 20. 

PACS numbers: 02.70.-c, 05.20.Gg, 82.20.Pm 



Path sampling has become an important tool to study 
rare events that are inaccessible for straightforward 
molecular dynamics (MD). Whereas the original path 
sampling approach [1] used a Monte Carlo sampling of 
trajectories with fixed lengths, the efSciency has been 
improved considerably by the introduction of the new 
transition interface sampling (TIS) technique that al- 
lows flexible path lengths. Besides a reduction in the 
required MD steps, it also yields a faster convergence by 
counting only effective crossing events. The TIS method 
has been applied to various systems ranging from pro- 
tein folding fsl] to nucleation [4|. Similar to the reactive 
flux (RF) approach , the TIS methods allow to deter- 
mine rate constants in terms of microscopic properties 
that do not sensitively depend of the choice of reaction 
coordinate (RC) and stable state definitions. However, 
whilst the efficiency of the RF methods drops dramati- 
cally whenever the RC does not capture the exact tran- 
sition mechanism, the TIS efficiency is relatively insensi- 
tive to the 'quality' of the RC 0]. This is an important 
advantage in high-dimensional complex systems where 
good RCs can be extremely difficult to find. TIS has also 
initiated the development of some new algorithms such as 
the partial path TIS (PPTIS) and forward flux sam- 
pling (FFS) . PPTIS uses a Markovian approximation 
to reduce the path length even further. FFS was espe- 
cially developed to deal with stochastic non-equilibrium 
systems. Similar to RF, the efhciency of these meth- 
ods is more sensitive to the RC Q. The advantageous 
scaling of TIS relies partly on the fact that it is an im- 
portance sampling on the dynamical factor. When this 
factor is low, direct evaluation as in RF becomes pro- 
hibitive. Secondly, due to the global character of trajec- 
tories, hysteresis effects are less likely to occur in path 
space than in phase space [6|. Finally, the non-locality 
of the shooting move might in principle allow the sam- 
pling of multiple reaction channels. However, if the reac- 
tion channels are very distinct, it might still take a long 
time for the shooting move to find these channels. In 
this letter, we introduce an additional technique based 
on replica-exchange methods to address this problem 
and test this method on the denaturation dynamics of 



the Peyrard-Bishop-Dauxois (PBD) model for DNA. 

The TIS algorithm works as follows. The first step is 
to define a RC and a set of related values Aq, Ai, . . . , A„ 
with Ai < Ai+i. The subsets of phase- or configuration 
points for which the RC is exactly equal to Ai basically 
define multidimensional surfaces or interfaces. These val- 
ues/interfaces should obey the following requirements: if 
the RC is lower than Ao = A/i, the system should be in 
the reactant state A; if the RC is higher than A„ = \b 
the system should be in the product state B; n and the 
positions for the interfaces in between should be set to op- 
timize the efficiency. Furthermore, the surface Xj^ should 
be set in such a way that whenever a MD simulation is 
released from within the reactant well, this surface is fre- 
quently crossed. The TIS rate expression can then be 
formulated as 

n-l 

kAB = /aPa(AbIAa) = /a n 'PA{X^+l\X^) (1) 

i=0 

Here, Ja is the flux through the first interface and can 
be computed by straight-forward MD. T^aIAbIAa) = 
'Pa(A„|Ao) is the probability that whenever the surface 
Aa is crossed. As will be crossed before Xa- The fac- 
torization of Va{Xm\Xo) into probabilities 7-'A(Ai+i |Ai) 
that are much higher than the overall crossing probabil- 
ity, is the basis of the importance sampling approach. It 
is important to note that VAiXi+i\Xi) are complicated 
history dependent conditional probabilities. If we con- 
sider all possible pathways that start at Xa and end by 
either crossing Xa or Xb , while having at least one cross- 
ing with Ai in between, the fraction that also crosses A^+i 
equals 7^A(Ai+i |Ai). This basically reduces the problem 
to a correct sampling of trajectories that should obey the 
Ai crossing condition. An effective method to achieve this 
is the so-called shooting algorithm [l[ . This Monte Carlo 
method randomly picks a time slice from the old existing 
path and makes a slight modification of this phase point. 
Then, this new phase point is used to propagate forward 
and backward in time yielding a new trajectory. In TIS, 
this propagation is stopped whenever the system enters A 
or B or, equivalently, whenever Ao or A„ are crossed. The 
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pathway is then accepted only if the backward trajectory 
ends in A and the total trajectory has at least one cross- 
ing with Xi. In addition, correct detailed balance rules 
are applied for the energy and path length fluctuations. 
The final result follows from the outcomes of a series of 
independent simulations {[md], [0+], . . . , The 
first is the MD simulation to compute The next ones 
are the path sampling simulations where indexes the 
surfaces that has to be crossed. In the same spirit of par- 
allel tempering (replica exchange) we could henceforth 
try to exchange the paths from one ensemble to the other, 
while running all these simulations simultaneously. The 
idea was already suggested in for PPTIS, but this 
is the first time that we show its effectiveness for TIS 
by making a small change to the algorithm. In order to 
have a full flexibility of swapping moves at all levels, we 
actually replace the MD simulation by another path sim- 
ulation [0~]. This ensemble consists of all possible paths 
that start at A^, then go initially in the negative direc- 
tion and end at the same interface A^. The flux /a can 
now be determined from 



(2) 



where (ip^th), (^path) ^'"'^ average path lengths in the 
[0^] and [0+] path ensembles respectively. The TIS al- 
gorithm is then as follows. At each step it is decided 
by an equal probability whether a series of shooting or 
swapping moves will be performed. In the first case, all 
simulations will be updated sequentially by one shooting 
move. In the second case, again an equal probability will 
decide whether the swaps [0~] ^ [0+], [1+] ^ [2+],... 
or the swaps [1+] <-> [2+], [3^] ^ [4+], . . . are performed. 
Each time that [0~] and [{n — 1)+] do not participate in 
the swapping move they are left unchanged. Also when 
the swapping move does not yield valid paths for both en- 
sembles, the move is rejected for the two simulations and 
the old paths are counted again. Note that the swap- 
ping moves do not require any force calculations. The 
only exception is [0~] ^ [0+]. Here, the last time step 
of the old path in the [0^] ensemble is used as initial 
point to generate a new trajectory in [0+] by integrating 
the equation of motion forward in time. Conversely, the 
initial point of the old path in [0+] is followed backward 
in time to generate a path in [0^]. The two types of 
swapping moves are illustrated in Fig. [TJ 

To test the efficiency of this new algorithm, we applied 
the TIS method with and without swapping moves to 
the denaturation transition of DNA usin g th e mesoscopic 
Peyrard-Bishop-Dauxois (PBD) model The PBD 

describes the DNA molecule as an one-dimensional chain 
of effective atom compounds yielding the relative base- 
pair separations j/i from the ground state positions. The 
total potential energy U for an N base-pair DNA chain is 
then given by C/(2/^) = Vi{yi)+j:f^^V,iy,)+Wiy,,y,.i) 
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FIG. 1: (color online) Illustration of the swapping move. The 
picture shows four possible pathways on a free energy sur- 
face corresponding to the [0~], [0^], [1^], [2^] ensembles. In 
the next step the swaps [0~] <-» [0^] and [l"""] « [2"*"] are 
performed simultaneously yielding four new pathways. Note 
that both [1"*"] and [2+] have moved to another reaction chan- 
nel. The alternative swapping move [0"*"] [1^] would have 
yielded a rejection as the [0^] does not cross Ai. 



with y^ = {yi} the set of relative base pair positions and 

1^.(2/0- A (3) 
IF(2/„ y,_i) = ^k(i + pe-^^y^+y-^A [y, - y,_i)2 



The first term Vi is the Morse potential describing the 
hydrogen bond interaction between bases on opposite 
strands. Di and ai determine the depth and width of 
this potential for the AT and GC base-pairs. The second 
term W is the stacking interaction. All interactions with 
the solvent and the ions are effectively included in the 
force-field. The constants K, p, a, Pat, DgCi QaTi Qgc 
were parameterized in Ref. [l^]. For finite chains, the 
associated state of the molecule is metastable as would 
be expected from a double stranded DNA molecule in a 
infinite solution. However, in order to fully dissociate, 
all base-pairs should reach the plateau of the Morse po- 
tential. As long as one base pair is still in the stack, 
it will likely pull all the others back to the associated 
state. At ambient conditions, the denaturation is a rare 
event and provides an excellent test for our methods. As 
X{y'^) = min[{?/i}] can describe both the associated and 
the disassociated state, we used this RC. For describing 
the reaction mechanism, a collective variable that de- 
pends on the positions of all particles might seem more 
appropriate. However finding the best possible RC is not 
the aim of this paper. Moreover, the choice of this RC 
has an additional advantage since it allows a very efficient 
third method that can be used as reference (see below). 

We consider a 20 AT base-pair DNA molecule that in- 
teracts with a 300 K Langevin thermostat with 7 = 50 
ps^^ to mimic aqueous damping. The time step was 
Ai = 1 fs and base pair masses were 300 amu. For the 
path sampling simulations we used aimless shooting [isj 



3 



where the velocities of the picked time shce are coni- 
pletely regenerated from MaxweUian distribution. This 
ensures a higher decorrelation between accepted paths 
than the standard shooting move where the velocities are 
only slightly changed. As the acceptance rates remained 
moderate > 0.3, the aimless shooting was found to be 
more efficient for this system. After a few trial simula- 
tions, all interfaces were positioned to have the optimum 
Pa(Aj+i|Aj) « 0.2 for all i 0, [IH. After this initializa- 
tion, all interface positions were fixed to Aq = 0, Ai = 
0.03, A2 = 0.07, A3 = 0.13, A4 = 0.21, A5 = 0.34, Xq = 0.7, 
and A7 = 1 A. In the next step, intensive calculations 
were run for each simulation which consisted of 10^ sim- 
ulation steps for the MD simulations and 4 • 10^ cycles 
for each path simulation. A technique to avoid complete 
separation in the MD simulation was applied ^l3| . Then, 
we repeated the same simulations with 2 • 10^ cycles using 
a 50 % swapping probability. The results are shown in 
tabled 



TABLE L Results of TIS simulations with and without path 
swapping. Errors are obtained using block averaging. 





Standard TIS 


Path Swapping 




value 


error (%) 


value 


error (%) 


/A(ns-i) 


304.9 


5.2 


291.8 


2.2 


Pa(Ai|Ao) 


0.256 


9.2 


0.249 


1.2 


Va(M^i) 


0.244 


4.6 


0.269 


1.9 


Pa{\3\\2) 


0.246 


5.0 


0.247 


2.4 


Va(\4\\3) 


0.310 


2.3 


0.307 


1.8 


Pa{\5\U) 


0.310 


1.3 


0.309 


1.2 


Pa (As 1 As) 


0.210 


1.1 


0.214 


1.3 


Pa(A7|A6) 


0.533 


0.5 


0.534 


0.8 


Pa(AsIAa) 


0.000165 


11.5 


0.000179 


4.2 


kAB (ns~^) 


0.0492 


12.6 


0.0524 


4.7 


2nd average 


0.0535 


22.4 


0.0533 


10.3 




FIG. 2: (color online) The matched overall crossing proba- 
bilities for the standard TIS and the parallel path swapping. 



Table H] shows that the parallel path swapping simula- 
tions have considerable lower errors despite the fewer sim- 
ulation cycles. Also, the construction of the overall cross- 
ing probabilities in Fig. [5] shows a much better matching. 
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FIG. 3: (color online) Error analysis of [0"*"] with and with- 
out swapping, a): the total running average "^^i^i/^ with 
Xi the j-th measurement, b) block average errors e{l) = 



E^(^j (0 - xy/x^M{M - 1) with Xj the average of the 
j-th block of length I, x the total average and M the num- 
ber of blocks. Horizontal lines at the plateaus indicate the 
actual relative error e. The total correlation is obtained by 



The final error was obtained by the error propagation rule 
Etot = \/Yl, ^i- the validity of these propagation rules 
is questionable for the correlated path swapping simula- 
tions, we also calculated the final results and errors in a 
different way. Here, we divided the simulations results 
in five blocks to obtain approximately independent rate 
constants. The errors, obtained by the standard devia- 
tion, are indeed higher. However, the errors for standard 
TIS increased as well by a similar factor. Probably, a 
much higher accuracy is needed to understand the effect 
of covariance terms for the errors in the swapping algo- 
rithm. We will come back to these results after we have 
compared the final outcome with that of a third method. 
This method is a very accurate implementation of the RF 
method that works due to some special characteristics of 
the system. Like in standard RF theory we write 

kAB = PA(AB)(A(2/^)x[path(y^)]){^(,«)=^^j (4) 

where Pa(Ab) is the probabihty that X{y^) — Xb given 
that the system is in state A. The second term is the (un- 
normalized) transmission coefficient and is calculated by 
releasing dynamical trajectories starting from the surface 
A B and x is a functional of the trajectory which corrects 
for fast recrossing events. Although different forms of x 
are used, the effective positive flux expression has shown 
to be the most efficient [^, ■ Here x is equal to 1 (0 
otherwise) only if A > and if the backward trajectory 
crosses Xa before Xb- The probabihty Pa (As) can be 
written as 



Pa (As) = 



/dy^^(A(y^)-AB)e-^^fa") 
J dy^9{XB - XiyN))e-0Uiv") 

j:jdy''S{X{y,)^XB)Uj^,d{y: 



(5) 



/dy~(l-nfe%A-As))e-/3^fe") 

As the integrals of Eq. (O are all of a special factorial 
form, we can apply the direct integration method of 14 1 



4 



yielding a precision of several digits. In the next step, we 
need to generate a representative set of configurations on 
the surface A^- It can be shown that ensemble averages 
with potential U{y^) and fixed min[{yi}] = As, are in 
fact identical to that of a freely moving chain on a trans- 
lational invariant potential U' that is related to U by 
U'{y^) = U{y^ — X{y^) + As). Hence, we can generate 
the required surface points by running a MD simulation 
using U'{y^), save every 1000th time step to dissolve 
correlations, and shift these configurations to the surface 
Xb ■ From these points, we release trajectories using nor- 
mal potential U and calculate x- We applied this method 
using a numerical integration step of dy = 0.01 A yield- 
ing the result Pa{Xb) = 5.316 • 10"^ A^^. Then we re- 
leased 4 • 10^ trajectories for which the initial points were 
generated using the dynamical shifted potential U'{y^). 
The transmission coefficient yielded 9.854 ± 0.066 A/ns 
and the combined result Uab = 0.0524 ns~^, which is in 
excellent agreement with the TIS results. 

Now we come back to the results of table H] and try to 
express the efficiency into the so-called efficiency times 
Tcff that are defined as the number of force calculations 
required to obtain a statistical error equal to 1. For the 
simulations [0+], [1^], . . . the efficiency times are given 
by i 

rij^l = I^^^UM^ (6) 

Pi 

Here, pi = VA{Xt+i\Xi) and Li = (ipa'^'h)/At which are in 
principle independent from the simulation method, is 
the ratio between the average cost of a simulation cycle 
and Li . Mi is the effective correlation. The results of the 



increase in efficiency of a factor 150. Inspection of Fig. [3]- 
a) reveals large fluctuation in the overall running average 
for standard TIS even after 4-10^ cycles. In contrast, the 
swapping results shows a much faster convergence. This 
is also reflected in the block-error analysis of Fig.[3]-b). 

The overall efficiency time is derived from t^s — e^Tsim- 
Here e the relative error in Uab and Tsim is the total sim- 
ulation time. The two ways of averaging show that the 
swapping moves give an overall improvement of approxi- 
mately 20. However, it is important to realize that apply- 
ing the same number of cycles for each simulation does 
not give the best possible performance. From the results 
of y, one can show that the efficiency can be improved 
by a factor of 7 if the optimal ratio of cycles proportional 

to oc \J l^iLi is applied. However, we believe that 
the path swapping efficiency can be improved by a simi- 
lar factor if we change the algorithm to allow an unequal 
distribution of shooting moves among the simulations. 
We are now working on such algorithms. 

To conclude, we have shown that TIS combined with 
path swapping can give a huge improvement of efficiency. 
For the denaturation of the PBD model of DNA, we ob- 
tained an improvement of approximately a factor 20. In- 
dividual path simulations were improved upto two orders 
of magnitude. Therefore, we believe that parallel path 
swapping can become an important method in any type 
of rare event simulations. 

I would like to thank Paolo Pescarmona for carefully 
reading this paper. 



TABLE II: Efficiency analysis 





Standard TIS 


Path Swapping 


L 


i 


TV 


Tcff (10^) 


L 




M 


roff(105) 


[md]/[0-] 


1 


1 


824 


27 


3262 


0.87 


146 


27 


[0+] 


115 


1.16 


11519 


45 


108 


0.83 


95 


0.3 


[1+1 


289 


1.05 


2087 


20 


325 


0.53 


261 


1 


[2+] 


764 


1.01 


3284 


78 


765 


0.51 


377 


4 


[3+] 


1832 


0.98 


921 


37 


1827 


0.49 


272 


5 


[4+] 


3768 


0.94 


327 


26 


3776 


0.47 


139 


6 


[5+] 


7464 


0.87 


121 


29 


7483 


0.43 


97 


11 


[6+] 


14391 


0.70 


109 


10 


14340 


0.35 


147 


6 


overall 
2nd aver. 








14907 
46660 








648 
2486 



efficiency analysis are given in table |TT] and show that 
the swapping moves decrease both ^ and TV. The effi- 
ciency times are all lowered by at least a factor of 2 for 
all path simulations. Spectacular is the decrease of corre- 
lation from 11519 to 95 in the [0+] simulation yielding an 



[1] C. Dellago, P. G. Bolhuis, and D. Chandler, J. Chem. 

Phys. 108, 9236 (1998). 
[2] T. S. van Erp, D. Moroni, and P. G. Bolhuis, J. Chem. 

Phys. 118, 7762 (2003). 
[3] P. G. Bolhuis, Proc. Nat. Acad. Sci. USA 100, 12129 

(2003). 

[4] D. Moroni, P. R. ten Wolde, and P. G. Bolhuis, Phys. 
Rev. Lett. 94, (2005). 

[5] D. Frenkel and B. Smit, Understanding molecular simu- 
lation, 2nd erf. (Academic Press, San Diego, CA, 2002). 

[6] T. S. van Erp, J. Chem. Phys. 125, (2006). 

[7] D. Moroni, P. G. Bolhuis, and T. S. van Erp, J. Chem. 
Phys. 120, 4055 (2004). 

[8] R. J. Allen, P. B. Warren, and P. R. ten Wolde, Phys. 
Rev. Lett. 94, 018104 (2005). 

[9] E. Marinari and G. Parisi, Europhysics Lett. 19, 451 
(1992). 

[10] T. Dauxois, M. Peyrard, and A. R. Bishop, Phys. Rev. 

E 47, R44 (1993). 
[11] T. S. van Erp and P. G. Bolhuis, J. Comput. Phys. 205, 

157 (2005). 

[12] A. Campa and A. Giansanti, Phys. Rev. E 58, 3585 
(1998). 

[13] B. Peters and B. L. Trout, J. Chem. Phys. 125, (2006). 
[14] T. S. van Erp, S. Cuesta-Lopez, J.-G. Hagmann, and M. 
Peyrard, Phys. Rev. Lett. 95, 218104 (2005). 



[15] J. B. Anderson, J. Chem. Phys. 62, 2446 (1975). 



